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[i]  Wc  show  that  ocean  temperature  profiles  can  be 
aeeurately  recovered  using  only  aeoustic  methods 
employed  at  the  sea  surface.  Using  a  towed  air  gun  array 
and  a  hydrophone  streamer,  thermohaline  boundaries  are 
ensonified  at  a  suite  of  frequencies  and  angles,  yielding 
travel  time  trajectories  and  reflectivities.  These  data  are 
inverted  via  full  waveform  inversion  to  estimate  sound 
speed  and,  subsequently,  a  temperature  profile.  The  high 
lateral  data  density  of  the  seismic  technique  offers  the 
potential  of  acoustically  derived  temperature  profiles  to  be 
used  to  constrain  models  of  oeean  mixing  and  internal 
waves.  Results  on  realistic  synthetic  data  show  that  sound 
speed  ean  be  recovered  with  arbitrary  aeeuraey  when  using 
broadband  data,  with  known  souree  function  and  recording 
geometry.  Application  to  field  seismic  data  (corroborated  by 
expendable  bathythermograph)  shows  that  even  with  a 
seismic  acquisition  system  not  specifically  calibrated  for 
seismic  oceanography,  temperature  contrasts  within  the 
ocean  can  be  recovered  to  within  one  degree  Celsius. 
Citation:  Wood.  W  T  ,  W.  S.  Holbrook,  M  K.  Sen,  and  P.  L. 
Stoffa  (2(X)8),  Full  waveform  inversion  of  refleciion  seismic  daia 
for  ocean  temperature  profiles,  Geophys.  Res.  Lett.,  35,  L04608, 
doi:  1 0. 1 029/2007GL032359. 

1.  Background  and  Motivation 

[2]  Marine  reflection  seismology  has  reecntly  been 
shown  capable  of  producing  detailed  images  of  oceanic 
thermohaline  finestructurc  at  lateral  resolution  of  ~10  m 
[Holbrook  et  al 2003].  Images  of  finestructure  have  been 
produced  in  numerous  ocean  settings,  including  fronts 
[Holbrook  et  al  ,  2003;  Noguchi  et  al.,  2006;  Tsuji  et  al., 
2005;  White  et  al .,  2006],  Meddies  [Klaeschen  et  al.,  2006; 
Pinheiro  et  al.,  2006],  intrathcrmocline  lenses  [ Bullock , 
2006],  warm-core  rings  [Seymour  et  al.,  2006],  water-mass 
boundaries  [Huthnance  et  al,  2006;  Nandi  et  al,  2004],  and 
thermohaline  staireascs  [Nandi  et  al ,  2006].  While  the 
images  alone  have  intrinsic  interest,  a  key  topic  of  eurrent 
research  in  “seismic  oceanography”  is  the  extraction  from 
seismie  data  of  quantitative  information  on  physical  ocean¬ 
ographic  processes  and  properties,  such  as  intcmal-wavc 
spectra  [Holbrook  and  Fer,  2005;  Krahmann  et  al,  2006] 
and  temperature  contrasts  [Paramo  and  Holbrook ,  2005]  In 
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this  paper  we  explore  the  possibility  of  inverting  seismic 
reflection  data  for  temperature-depth  profiles  in  the  ocean, 
using  one-dimensional  (ID)  full-waveform  inversion,  dem¬ 
onstrated  on  seismic  data  eo-loeated  with  XBT  profiles  in 
the  Norwegian  Sea. 

[3]  Seismie  oceanography  data  arc  especially  well  suited 
to  one -dimensional  waveform  inversion  approaches,  be¬ 
cause  thermohaline  boundaries  in  the  ocean  are  very  nearly 
flat  and  horizontal,  lateral  variations  in  sound  speed  arc 
small,  no  converted  shear  waves  are  present,  and  interbed 
multiples  are  negligible  due  to  the  small  reflection  coef¬ 
ficients  (^0  001).  Through  full  waveform  inversion  [e.g,, 
Singh  et  al,  1993;  Wood  et  al,  1994;  Korenaga  et  al, 
1997],  every  reflection  within  the  water  column  is  mod¬ 
eled  simultaneously,  resulting  in  a  1-D  profile  of  sound 
speed,  which  can  then  be  readily  converted  to  temperature 
via  an  equation  of  state  [e  g.,  Chen  and  Millcro,  1977], 
Although  ihc  inversion  algorithm  is  capable  of  distin¬ 
guishing  independently  varying  density  and  sound  speed, 
density  contrasts  are  typically  much  smaller  than  sound 
speed  contrasts  in  the  ocean  and  contribute  little  to  the 
refleetanee  in  the  Norwegian  Sea  data  set  used  here 
[Paramo  and  Holbrook,  2005],  so  we  assume  for  this 
study  that  all  reflectivity  in  the  water  eolumn  is  associ¬ 
ated  with  sound  speed  contrasts. 

2.  Inversion  Method 

[4]  The  groundwork  for  the  class  of  problems  and  sol¬ 
utions  used  here  has  been  diseussed  extensively  by  Menke 
[1989]  and  Tarantola  [1987],  and  the  application  to  reflec¬ 
tion  seismic  data  is  based  on  the  work  of  McAulay  [1985, 
1986],  Dietrich  and  Kormendi  [1990],  Amundsen  and  Ursin 
[1991],  and  Wood  et  al.  [1994].  We  define  the  general 
forward  problem  as 

d  =  g(m)+v  (I) 

where  m  is  a  veetor  of  nmod  model  parameters  (in  this  ease 
sound  speeds),  g  is  a  non-linear  operator,  in  this  case 
containing  information  on  acoustic  wave  propagation  and 
the  experimental  configuration,  v  is  a  veetor  of  additive 
noise,  and  d  is  a  vector  of  ndal  data  (in  this  case  waveform 
amplitudes)  that  would  be  observed  if  the  water  column 
sound  speed  was  perfectly  described  by  m,  if  g  was  a 
perfect  theoretical  relation,  and  if  v  was  identically  zero. 
(Vectors  and  vector  functions  of  vectors  are  cast  in  lower 
case  bold  type  while  matrices  are  east  in  uppercase  bold 
type)  The  inverse  problem  ean  then  be  defined  as 

m  =  y(d  —  v),  (2) 
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where  v  ean  be  absorbed  into  the  data  vector,  and  y  is  a 
generalized  inverse,  the  structure  of  which  is  the  subjeet  of 
much  discussion  in  the  literature,  (e.g.,  monographs  by 
Tarantola  [1987]  and  Menke  [1989]).  For  linear  problems 
and  wisely  ehosen  objective  functions  a  single  evaluation  of 
the  slope  and  curvature  of  the  objective  function  at  any 
point  can  be  used  to  generate  the  best  possible  estimate  of 
the  location  of  the  objective  function  minimum  and  optimal 
model.  If  the  problem  is  almost  linear,  then  an  evaluation  of 
the  slope  and  curvature  result  in  an  improved,  but  not 
optimal  model  estimate,  and  may  have  to  be  re-evaluated, 
making  the  solution  iterative,  but  still  linearizable  (as  is  the 
case  for  sound  speed  in  this  study).  We  use  the  least  squares 
objective  function,  S,  defined  in  vector  notation  as; 

2S[Ad'C/)1  Ad  +  Am'Cj  Ain], 

with 

Ad  =  [d,Vfl  -  dohs]  -  [g(m)„  d«*T], 

and 

Am  =  [m(J  -  m0]  (3) 

where  CD  and  Cm  are  prior  data  and  model  eovarianee 
matrices  respectively,  and  m0  and  mn  are  model  vectors 
containing  the  a  priori  and  nth  trial  model  parameters, 
respectively,  (superscript  t  denotes  transpose).  The  veetor 
dohs  contains  the  observed  data,  and  dsyn  contains  the 
synthetie  data  generated  from  the  n  h  model  iteration  by 
the  operator  g.  The  vector  dsyn  ean  also  be  expressed  as  the 
vector  function,  g(mn).  The  objective  function  consists  of 
both  a  model  and  data  error,  and  it  is  the  combination  of 
these  errors  that  is  to  be  minimized.  Even  if  the  data  are 
matehed  exactly  the  objective  function  may  be  quite  large  if 
the  corresponding  model  is  far  from  the  prior  model.  Errors 
in  the  theory  of  the  forward  modeling  g(mn)  take  the  same 
form  as,  and  are  absorbed  into,  the  data  covariances  in  CD. 
For  this  study  these  modeling  errors  oecur  when  the  ocean 
structure  is  not  1-D  (dipping  or  discontinuous  layers),  not 
isotropic  (wave  speeds  are  directionally  dependent),  or 
inelastic  (wave  amplitudes  are  attenuated  as  a  function  of 
time). 

[5]  For  this  study  we  chose  the  least  squares  error 
solution  defined  above  with  Newton’s  method  of  minimi¬ 
zation.  The  convergence  for  linear  problems  is  quiek,  the 
mathematical  foundation  for  the  solution  is  well  known 
[ Tarantola ,  1987;  Menke ,  1 989],  and  a  quantitative  estimate 
of  the  variance  (uncertainty)  in  the  solution  can  be  easily 
obtained  from  the  curvature  of  the  error  surface  at  the  point 
of  minimum  error. 

[6]  The  iterative  form  of  New  ton's  method  can  be  written 
in  multi-dimensional  form  [ Tarantola ,  1987]  as 

m„,i  m„  -  [(^S/<9m:]„_l  (i9S/9m]n  (4) 

where  d  denotes  the  partial  derivative,  and  S  is  the  value  of 
the  objective  function.  Here  both  the  slope  and  curvature  of 
the  model  space  at  the  nth  iteration  are  used  to  find  the  next, 
or  n  +  1th  model  estimate. 


[7]  Taking  the  derivative  of  S  to  find  the  local  multidi¬ 
mensional  slope  yields 

[aS/am]B  =  GJ.C0  Ad  +C'1  Am  (5) 

where  the  matrix  Gn  (sensitivity  or  Frechct  derivative 
matrix)  contains  the  sensitivity  of  eaeh  data  parameter  to 
each  model  parameter 

Gt  =  [<W7cW]„  (6) 

[a]  Taking  the  seeond  derivative  to  find  the  local  multi¬ 

dimensional  curvature  yields; 

rs  G'C^'G,  +  C"1 ,  (7) 

where  the  neglected  term  is  small  when  the  forward 
problem  is  nearly  linear 

[9]  Combining  equations  (4),  (5),  and  (7)  gives  the 
expression  used  in  this  study, 

"Wi  =  m„  +  [G'nCD[G„+Cm'}  1  [G'C,,1  Ad  +  C,,1  Am],  (8) 

where  Gn  is  re-computed  at  each  iteration,  and  each  of  the 
other  components  on  the  right  hand  side  is  known  at  the 
start  of  the  inversion.  The  nmod  x  nmod  matrix,  [G„  C^1  Gn  + 
Cm1],  referred  to  as  the  Hessian,  is  solved  via  singular  value 
decomposition  (SVD)  [e.g.,  Menke ,  1989]. 

3.  Application  to  Synthetic  Seismic  Data 

[10]  To  apply  equation  (8)  to  seismic  reflection  data,  we 
parameterize  the  problem  such  that  the  model  is  a  scries  of 
sound  speeds  corresponding  to  layers  in  the  water  column 
whose  time  thicknesses  are  held  constant  at  the  seismic  data 
sample  increment  of  0.004  seconds.  This  ensures  that  any 
layer  that  can  affect  the  seismic  data  can  be  completely 
modeled,  i.c.  any  synthetic  seismogram  can  be  reproduced 
exactly.  We  parameterize  the  data  as  frequency  domain, 
plane-wave  seismograms  as  obtained  by  a  Founer-Hankcl 
transform  of  a  common  midpoint  (CMP)  gather.  This  allows 
for  use  of  the  very  efficient  reflectivity  method  \Kennett  and 
Kerry ,  1979]  as  the  forward  operator  g(m)  in  equation  (1). 
Our  use  of  the  reflectivity  method  also  requires  knowledge 
of  the  souree  function  (wavelet),  whose  shape  (phase)  was 
determined  by  iteratively  modeling  frequency  component 
phases  to  minimize  the  total  energy  of  a  senes  of  field  data 
traces  [Wood,  1999], 

[n]  The  degree  to  which  each  of  the  data  and  model 
misfits  are  minimized  is  controlled  by  the  a  priori  cova¬ 
riances  in  Cm  and  CD.  For  the  trials  presented  in  this  study 
both  Cm  and  Cq  are  assumed  to  be  purely  diagonal 
matrices,  with  eaeh  diagonal  a  constant  value.  The  data 
and  model  covariances,  which  correspond  to  the  squared 
uncertainties,  were  ehosen  as  0.01  and  300  respectively. 
Note  that  these  parameters  can  also  be  regarded  as  regular¬ 
ization  parameters  that  introduce  stability  in  the  inversion. 
The  chosen  values  result  in  the  data  misfit  being  much  more 
heavily  weighted  than  the  model  misfit,  appropriate  in  this 
ease  where  we  are  mueh  less  certain  about  the  a  priori 
model  than  about  the  observed  data. 
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[12]  To  test  the  sensitivity  of  the  inversion  we  generate  a 
synthetic  seismogram  based  on  the  same  Norwegian  Sea 
seismic  data  and  coincident  XBT  (expendable  bathyther¬ 
mograph)  profiles  presented  by  Paramo  and  Holbrook 
[2005],  displayed  here  in  the  intercept  time-slowness 
(Tau-p)  domain,  rather  than  the  frequency-slowness  (u>p) 
domain  for  ease  of  interpretation  (Figure  lb).  The  temper¬ 
atures  from  the  XBT  were  converted  to  sound  speed 
assuming  a  constant  salinity  of  32  PSU,  and  using  these 
sound  speeds  were  converted  to  time  and  re-sampled  at  4 
millisecond  intervals  (Figure  la).  Errors  in  overall  salinity 
manifest  only  at  the  lowest  frequencies  (less  than  1  Hz),  and 
because  an  a  priori  starting  model  supplies  this  frequency 
band,  the  sensitivity  of  the  inversion  to  overall  salinity  is 
small.  The  XBT  and  the  wavelet  estimated  from  the  data 
were  used  to  generate  the  synthetic  seismogram  shown  in 
Figure  lb. 

[13]  Two  inversion  results  of  ideal  synthetic  data  arc 
presented  in  Figures  la-le,  using  starting  models  (red 
curves)  that  arc  0.0  and  5.0  Hz  low  pass  (cosine  tapered 
from  0  to  5  Hz)  versions  of  the  true  model,  and  both 
performed  over  the  entire  0-125  Hz  frequency  band  of  the 
seismic  data.  When  using  the  5  Hz  starting  model,  the 
algorithm  converged  to  the  true  model  (black  curve, 
Figure  la),  matching  the  data  to  within  one  part  in  a 
thousand,  leaving  only  a  low  amplitude,  low  frequency 
residual  (Figure  le).  Using  the  0  Hz,  constant  sound  speed 
of  1499  m/s  as  the  starting  model,  the  algorithm  accurately 
recovered  the  higher  frequency  portions  of  the  model,  but 
recovery  of  the  lower  frequency  components  is  incomplete 
below  the  strongest  reflection  event  at  about  0.6  seconds 
(Figure  la),  The  degradation  is  due  mainly  to  the  laek  of 
strong  reflections  in  this  portion  of  the  data  set,  whose 
trajectories  drive  the  recovery  of  the  low  frequency  model 
components,  henee  the  low  frequency  residual  in  Figure  lc. 
The  higher  frequency  components  of  the  profile  (detail) 
actually  enable  the  recovery  of  the  lower  frequency  com¬ 
ponents  (background)  through  the  generation  of  reflection 
events. 

[14]  For  the  ‘  noiseless”  synthetie  examples  in  Figures 
la-le,  the  aeeuraey  of  the  inversion  is  limited  only  in  that 
the  ratio  of  the  smallest  to  the  largest  eigenvalue  in  the 
Hessian  matrix  must  be  larger  than  the  machine  precision. 
This  places  the  effective  noise  floor  at  the  level  of  maehme 
precision,  i.c.  10  K,  for  all  synthetic  trials  in  this  study, 
allowing  exceptionally  low  amplitude  frequency  compo¬ 
nents  to  contribute  to  recovering  the  thermal  profile.  Re¬ 
cently  reported  techniques  may  eliminate  the  need  to 
compute  the  matrix  inverse,  allowing  even  greater  accuracy 
[Sen  and  Roy ,  2003;  Roy  et  al.t  2005]. 

[15]  The  successful  recovery  of  the  correct  temperature- 
depth  profile  by  the  seismic  waveform  inversion  is  not 
merely  confirmation  that  the  inversion  technique  and  algo¬ 
rithm  work,  but  rather  a  demonstration  that,  for  seismie  data 
with  typical  frequency  range  of  10-125  Hz,  realistic  fines- 
trueture  ean  be  fully  resolved  to  the  extent  that  eertain 
practical  conditions,  whieh  we  discuss  next,  arc  met. 

4.  Application  to  Field  Data 

[16]  We  next  apply  the  inversion  to  field  data  coincident 
with  an  XBT  so  the  inversion  technique  can  be  indepen¬ 


dently  corroborated.  The  data  were  originally  sampled  at 
0.002  milliseconds,  however,  most  of  the  source  band  was 
below  125  Hz,  so  we  low  pass  filtered  and  re-sampled  the 
data  to  0.004  milliseconds  This  reduces  the  number  of 
model  parameters,  nmod,  by  a  factor  of  two,  greatly  facili¬ 
tating  the  SVD  of  the  Hessian  matrix. 

[17]  Application  of  the  inversion  to  field  data  also  requires 
transforming  the  data  into  the  plane  wave  (u>p)  domain  [c.g., 
Bry’sk  and  McCowan ,  1986a,  1986b]  and  using  tapering  to 
minimize  the  transform  artifacts.  The  inversion  is  performed 
over  the  slowness  range  p  =  0-0.6  s/km,  or  approximately 
0-64  degrees  ineidenee  angle.  The  small  group  size  in  the 
array  affects  the  directivity  only  slightly  (1.5%  at  90  Hz  and 
65  degrees  incident  angle)  so  we  assume  it  is  negligible. 

[is]  Both  the  absolute  value,  and  the  angle  dependence  of 
the  reflection  coefficient  are  required  for  the  inversion.  To 
convert  seismic  signals  to  reflection  coefficients  we  com¬ 
pare  the  water  bottom  primary  and  first  sea  surfaec  multiple 
[1 Warner ,  1990]  resulting  in  a  scale  faetor  of  5.0  x  104  +/— 
2.0  x  104.  The  large  uncertainty  was  due  to  significant 
interference  between  the  multiple  and  primary  reflections 
from  subsurface  stratigraphy.  The  data  scaling  error  mani¬ 
fests  primarily  as  a  constant  multiplicative  factor  to  the 
entire  a  posteriori  model,  but  beeause  the  lowest  frequen¬ 
cies  in  the  a  posteriori  model  rely  almost  completely  on  the 
lowest  frequencies  in  the  a  priori  model,  the  inversion  at 
these  frequencies  is  mostly  insensitive  to  scaling  error. 
Within  the  data  band,  the  scaling  errors  manifest  as  an  error 
in  the  deviation  from  the  smooth  background,  which  can  be 
incorporated  into  the  model  uncertainties  associated  with 
the  smooth  starting  model.  Although  not  done  so  here,  the 
scale  faetor  eould  potentially  be  included  as  a  model 
parameter  for  which  to  invert. 

[19]  The  observed  data  in  the  plane  wave  domain  arc 
shown  in  Figure  le,  along  with  the  best  (smallest)  residual 
from  the  inversion  (Figure  If),  and  the  model  responsible 
for  the  best  residual  (Figure  Id).  As  in  the  synthetic  data 
example,  two  inversion  results  are  presented;  corresponding 
to  starting  models  that  are  low-pass  filtered  (eosine  taper 
from  0-30  Hz)  versions  of  the  XBT  profile  The  zero  Hz  case 
is  extreme  and  used  here  for  illustrative  purposes.  In  most 
eases,  analyzing  the  trajectory  of  the  sea  bottom  reflection 
can  yield  an  average  velocity  profile  in  the  0-1.0  Hz  range. 
The  30  Hz  starting  model  supplies  the  components  that  are 
poorly  constrained  or  missing  in  the  data  band  (about  25- 
80  Hz)  and  the  inversion  supplies  the  rest,  resulting  in  a 
low  data  residual  (Figure  If).  This  result  also  provides  a 
gross  eheek  on  the  inversion,  confirming  that  the  algorithm 
will  converge  on  what  we  expect  is  the  true  model,  (i.e.  the 
XBT  profile  lies  very  near  the  position  of  the  global 
minimum). 

[20]  Even  when  a  constant  sound  speed  is  used  as  the 
starting  model,  the  lowest  frequency  components  of  the 
model  are  well  recovered  down  to  the  major  reflection  at 
0.6  seeonds.  There  are  two  exceptions:  1)  the  shallow 
portion  corresponding  to  the  warm  surfaee  waters  where 
the  corresponding  data  were  muted  due  to  interference  with 
the  direet  wave,  and  2)  the  event  at  about  0.75  s  that  results 
from  a  more  gradual  thermal  transition  than  the  event  at 
0.6  s,  placing  it  near  the  low  end  of  the  spectrum  where 
recovery  is  weak.  The  inset  in  Figure  le  shows  that  the 
magnitude  of  the  temperature  step  at  0.6  s  is  recovered 
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Figure  1.  The  sound  speed  profile  from  (a)  an  XBT  cast  (true  model,  black  curve)  was  used  to  generate  (b)  a  synthetic 
data  set,  used  as  the  observed  data  in  a  test  of  the  inversion.  The  starting  models  (red  curves)  are  0  Hz  and  5  Hz  low  pass 
filtered  versions  of  the  true  model.  The  green  curve  is  the  inversion  from  the  0  Hz  starting  model,  and  the  inversion  result 
from  the  5  Hz  starting  model  (not  shown)  is  effectively  coincident  with  the  true  model,  resulting  in  (c)  a  data  residual 
(difference  between  observed  and  synthetic,  magnified  one  thousand  times)  that  is  extremely  small,  (e)  The  field  data  were 
acquired  coincidcntly  with  (d)  the  XBT  (black  curve),  which  can  be  used  to  independently  corroborate  the  inversion  result. 
The  starting  models  (sec  Figure  Id,  red  curves)  are  0  Hz  and  30  Hz  low  pass  filtered  versions  of  the  XBT  profile.  The  two 
inversion  results  arc  shown  in  green.  When  the  low  frequency  portion  of  the  model  is  supplied  in  the  starting  model,  the 
inversion  result  matches  the  XBT  and  the  (t)  data  residual  is  minimized. 
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Figure  2.  The  discrepancy  between  the  spectra  of  the  XBT 
and  inverted  sound  speed  profile  using  the  0  Hz  (blaek)  and 
30  Hz  (gray)  starting  model  shows  that  either  starting  model 
will  result  in  roughly  equal  recovery  of  frequencies  above 
about  30-40  Hz.  Supplying  frequencies  below  30  Hz  in  the 
sound  speed  model  (gray)  significantly  reduees  the 
discrepancy. 

accurately  regardless  of  the  starting  model,  but  recovery  of 
the  absolute  temperature  requires  aeeurate  background  (low 
frequency)  components.  Figure  2  shows  more  quantitatively 
how  spectra  of  the  inverted  sound  speed  profiles  differ  from 
the  spectrum  of  the  XBT.  For  the  profile  recovered  using  the 
30  Hz  starting  model,  the  spectral  difference  (gray  eurve)  is 
much  smaller  at  the  lower  frequencies,  because  these 
frequencies  have  been  supplied  by  the  starting  model.  The 
discrepancy  is  much  larger  (blaek  curve)  for  the  profile 
recovered  using  the  zero  Hz  starting  model.  Frequency 
components  well  within  the  data  band  (25-80  Hz)  are 
recovered  roughly  equally  well  regardless  of  the  starting 
model. 

[21]  Based  on  the  successful  recovery  of  the  low  frequen¬ 
cy  components  above  0.6  seconds,  and  the  results  of  the 
trials  on  synthetic  data,  we  would  expect  that  had  the 
seafloor  reflection  been  included  in  the  inversion,  the 
seismie  data  alone  would  have  had  the  low  frequency 
components  necessary  to  recover  the  absolute  temperature. 
These  components  would  not  need  to  have  been  contributed 
via  the  starting  model.  The  seafloor  reflection  was  not 
included  in  the  inversion  beeause  the  vast  difference  in 
amplitude  (about  3  orders  of  magnitude)  between  seafloor 
and  water  column  reflections  caused  significant  artifacts  in 
the  frequency  domain  processes  (filtering,  plane-wave  trans¬ 
form,  and  inversion)  used  here. 

[22]  Unlike  Ihe  synthetic  data  example,  there  are  small 
discrepancies  between  the  best-recovered  model  and  the 
“taie”  XBT  profile.  The  XBT  may  not  be  located  at  exactly 
the  same  position  as  the  CMP.  Also,  the  CMP  is  composed 
of  several  subsequent  source  firings,  acquired  up  to  several 
minutes  and  tens  of  meters  apart,  over  which  some  averag¬ 
ing  is  inevitable.  Further,  each  souree  fire  may  be  different 
due  to  irregularities  in  firing  time,  air  pressure,  sea  state,  or 
other  meehanieal  phenomena.  The  irregularities  are  usually 


negligible,  but  may,  along  with  any  artifaets  of  the  finite 
integral  transform,  manifest  either  as  errors  in  the  wavelet, 
which  ean  result  in  unwanted  deviations  in  the  final  model, 
or  as  artifaets  in  the  data  that  cannot  be  modeled  but  may 
draw  the  algorithm  away  from  its  desired  eourse  of  conver¬ 
gence.  Acquisition  irregularities  or  souree  wavelet  errors  are 
the  most  likely  eause  of  the  approximately  10  Hz  undu¬ 
lations  in  the  recovered  profiles  in  Figure  Id,  and  eould 
likely  be  significantly  mitigated  in  an  experiment  designed 
specifically  for  seismie  oceanography. 

5,  Conclusions 

[23]  We  have  demonstrated  that  full  waveform  inversion 
ean  recover  ocean  temperature  profiles  from  surface  tow'cd 
seismie  measurements  alone.  Inversions  of  synthetic  seis¬ 
mic  data  show  the  technique  ean  resolve  oceanic  finestrue- 
ture  at  the  5  m  vertical  scale  with  arbitrary  aeeuraey.  In  a 
test  on  field  data,  the  aeeuraey  of  the  recovered  profile 
depends  on  how  well  the  frequency  bands  of  the  seismie 
data  and  starting  model  eover  the  band  of  the  desired 
profile,  and  on  the  presence  of  reflections  throughout  the 
profile  to  provide  the  low  frequency  components  of  the 
profile.  The  latter  requirement  ean  be  relaxed  if  sufficient 
XBT  data  are  present  regionally  to  provide  background 
temperature-depth  profiles  for  starting  models.  The  recovery 
of  sharp  thermal  ehanges  is  generally  more  aeeurate  than  the 
recovery  of  the  absolute  temperature.  With  customary 
seismic  equipment  and  recording  techniques,  we  have 
estimated  the  magnitude  of  a  temperature  contrast  to  within 
approximately  20%,  or  0.5  degrees  C  of  its  value  as 
measured  directly  with  an  XBT.  This  aeeuraey  eould  be 
improved  significantly  w  ith  more  accurate  measurements  of 
the  souree  wavelet,  and  broader  band  data.  Our  rcsulls 
suggest  that  augmenting  sparse  direct  temperature  measure¬ 
ments  (e.g.,  XBTs)  with  inversions  of  seismie  data,  may 
result  in  an  effective  means  of  achieving  high  lateral  and 
vertical  resolution  thermal  eross-seetions  over  extensive 
regions  of  the  oeean. 

[24]  Acknowledgments.  The  authors  thank  the  reviewers  for  their 
helpful  comments.  This  work  was  funded  by  NSF-OCF  0452744  and  NSF 
OCE  0648620. 

References 

Amundsen,  L.,  and  B.  Ursin  (1991),  Frequency -wavenumber  inversion  of 
acoustic  data.  Geophysics,  56.  1027  1039. 

Brysk,  H.,  and  D  W.  McCowan  (1986),  A  slanl  stack  procedure  for  poinl- 
source  data.  Geophysics ,  51,  1370  1386. 

Brysk,  H.,  and  D.  W.  McCowan  (1986),  Edge  effects  in  cylindrical  slant 
stacks,  Geophys.  J  R.  Astron.  Soc .,  #7,  801  813. 

Bullock.  A.  D.  (2006),  Enhanced  fincstructurc  around  an  intrathcrmoclinc 
lens  in  the  Norwegian  Sea,  Eos  Trans.  AGO,  #7(36),  Ocean  Sci  Meet. 
Suppl .  Abstract  OS1 31-04. 

Chen,  C,  T  ,  and  F  J  Millcro  (1977),  Speed  of  sound  in  seawater  al  high 
pressures,  J.  Acoust.  Soc.  Am.,  62,  1 129-  1135. 

Dietrich,  M.,  and  F.  Kormcndi  (1990).  Perturbation  of  the  plane-wave 
reflectivity  of  a  depth-dependent  clastic  medium  by  weak  inhomogene  - 
ilies.  Geophys.  J.  Ini.,  100 ,  203  214. 

Holbrook,  W.  S.,  and  I.  Fcr  (2005),  Ocean  internal  wave  spectra  inferred 
from  seismic  reflection  lransccls,  Geophys.  Res.  Lett.,  32.  LI 5604, 
dot:  1 0. 1 029/2005 G  L0 2 3  7 3  3 . 

Holbrook,  W.  S.,  ct  al.  (2003),  Thcmiohalinc  fine  structure  in  an  oceano¬ 
graphic  front  from  seismic  reflection  profiling.  Science,  301,  821  824. 
Huthnancc.  J.,  cl  al.  (2006),  Imaging  the  time  dependence  of  water  bound¬ 
aries  in  the  Faroc-Shctland  Channel,  Eos  Tram.  AGV,  #7(36),  Ocean  Sci. 
Meet.  Suppl.,  Abstract  OS  131-06. 


5  of  6 


L 04608 


WOODS  ET  AL.:  WAVEFORM  INVERSION  OF  REFLECTION  SEISMIC  DATA 


1,04608 


Kcnnct,  B.  L.  N.,  and  N.  J.  Kerry  (1979),  Seismic  waves  in  a  stratified  half- 
space,  Geophys.  J.  R.  Astrvn.  Soc .,  57,  557  583. 

Klacschcn,  D.,  ct  al.  (2006),  Seismic  images  and  properties  of  a  Meddy,  Eos 
Trans.  AGV,  #7( 36),  Ocean  Sci.  Meet.  Suppl.,  Abstract  OS1 31-03. 

Korenaga,  J.,  W.  S.  Holbrook,  S.  C.  Singh,  and  T.  A.  Minshull  (1997), 
Natural  gas  hydrates  on  the  southeast  U.S.  margin:  Constraints  from  full 
waveform  and  travel  time  inversions  of  wide-angle  seismic  data,  J.  Geo¬ 
phys.  Res.,  102,  15,345  15,365. 

Krahmann,  G.,  ct  al.  (2006),  Internal  wave  energy  estimated  from  seismic 
reflection  data,  Eos  Trans.  AGV,  #7(36),  Ocean  Sci.  Meet.  Suppl.,  Ab¬ 
stract  OS  1 61-03. 

McAulay,  A.  D.  (1985),  Pre-stack  inversion  with  plane-layer  point  source 
modeling.  Geophysics,  50,  77  89. 

McAulay,  A.  D.  (1986),  Plane  layer  prestack  inversion  in  the  presence  of 
surface  reverberation.  Geophysics ,  51,  1789-1800. 

Menke,  W  (1989),  Geophysical  Data  Analysis:  Discrete  Inverse  Theory, 
Accdcmic,  San  Diego,  Calif 

Nandi,  P,  W.  S.  Holbrook,  S.  Pcarsc,  P.  Paramo,  and  R.  W.  Schmitt  (2004), 
Seismic  reflection  imaging  of  water  mass  boundaries  in  the  Norwegian 
Sea.  Geophys.  Res.  Lett.,  31,  L23311,  doi:10.l029/2004GL021325. 

Nandi,  P.,  ct  al.  (2006),  Seismic  imaging  of  a  thermohalinc  staircase  in  the 
western  tropical  Atlantic,  Eos  Trans.  AGIJ,  #7(36),  Ocean  Sci.  Meet. 
Suppl.,  Abstract  OS1 61-01 

Noguchi,  T.,  ct  al.  (2006),  Thermohalinc  interleaving  in  the  Kuroshio  Ex¬ 
tension  Front:  Seismic  reflection  imagery  and  in  situ  measurements,  Eos 
Trans.  AGV,  #7(36),  Ocean  Sci.  Meet  Suppl.,  Abstract  OSI 61-06. 

Paramo,  P.,  and  W  S.  Holbrook  (2005),  Temperature  contrasts  in  the  water 
column  inferred  from  amplitudc-vcrsus-offsct  analysis  of  acoustic  reflec¬ 
tions.  Geophys.  Res.  Lett.,  32,  L24611,  doi:l0.1029/2005GL024533. 

Pinheiro,  L.  M.,  ct  al.  (2006),  Detailed  2-D  imaging  of  the  Mediterranean 
Outflow  and  Meddies  off  W  Iberia  from  multichannel  seismic  data,  Eos 
Trans.  AGV,  #7(36).  Ocean  Sci.  Meet.  Suppl.,  Abstract  OSI 31-02. 

Roy,  L.,  M.  K.  Sen,  D.  Blankenship,  P.  L.  StofTa,  and  T.  Richter  (2005), 
Inversion  and  uncertainty  estimation  of  gravity  data  using  simulated 
annealing:  An  application  over  Lake  Vostok,  East  Antarctica,  Geophy¬ 
sics,  70,  JI-J12. 


Sen,  M.  K.,  and  1.  G.  Roy  (2003).  Computation  of  differential  seismograms 
and  iteration  adaptive  regularization  in  prestack  waveform  inversion. 
Geophysics,  68,  2026  2039. 

Seymour,  J.  C.,  ct  al.  (2006),  Velocity  structure  of  a  Loop  Current  eddy 
using  the  gradicnt-cun  ent  method  on  a  seismic  section  from  the  Gulf  of 
Mexico,  Eos  Trans.  AGV,  #7(36),  Ocean  Sci.  Meet.  Suppl.,  Ahstract 
OSI  31-05 

Singh.  S.  C.,  T.  A.  Minshull,  and  G  D.  Spence  (1993),  Velocity  structure  of 
a  gas  hydrate  reflector.  Science,  260,  204  207. 

Tarantola,  A.  (1987),  Inverse  Problem  Theory:  Methods  for  Data  Fitting 
and  Model  Parameter  Estimation,  Elsevier,  New  York 

Tsuji,  T.,  T,  Noguchi,  H  Niino,  T.  Matsuoka,  Y.  Nakamura,  H  Tokuyama, 
S.  Kuramoto,  and  N.  Bangs  (2005).  Two-dimensional  mapping  of  fine 
structures  in  the  Kuroshio  Current  using  seismic  reflection  data,  Geophys. 
Res.  Lett. ,  32,  L14609,  doi:!0.1029/2005GL023095. 

Warner,  M.  (1990),  Absolute  reflection  coefficients  from  deep  seismic  re¬ 
flections,  Techtonophvsics,  173 ,  15  23. 

White,  N.,  ct  al.  (2006),  Seismic  imaging  of  fronts  associated  with  Antarctic 
Circumpolar  Current  near  Drake  Passage,  Eos  Trans.  AGV,  #7(36), 
Ocean  Sci.  Meet.  Suppl.,  Abstract  OS  161 -04. 

Wood,  W.  T.  (1999),  Simultaneous  deconvolution  and  wavelet  inversion  as 
a  global  optimization.  Geophysics,  64,  1 108  1115. 

Wood,  W.  T.,  P.  L.  StofTa,  and  T.  H.  Shipley  (1994).  Quantitative  detection 
of  methane  hydrate  through  high-resolution  seismic  velocity  analysis. 
J.  Geophys.  Res.,  99,  9681  9695. 


W.  S.  Holbrook,  Department  of  Geology  and  Geophysics,  University  of 
Wyoming,  1000  East  University  Avenue,  Dept.  3006,  Laramie,  WY  8207 1  - 
3006,  USA. 

M  K.  Sen  and  P.  L.  Stoffa.  Institute  for  Geophysics,  University  of  Texas 
at  Austin,  Austin,  TX,  USA. 

W.  T.  Wood,  Code  7432,  Naval  Research  Laboratory,  Stenms  Space 
Center,  MS  39529,  USA.  (warren .wood(a;nrlssc  navy.mil) 


6  of  6 


